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Abstract 

We consider a model of large regulatory gene expression networks where the thresholds activating the 
sigmoidal interactions between genes and the signs of these interactions are shuffled randomly. Such an 
approach allows for a qualitative understanding of network dynamics in a lack of empirical data concerning 
the large genomes of living organisms. 

Local dynamics of network nodes exhibits the multistationarity and oscillations and depends crucially 
upon the global topology of a "maximal" graph (comprising of all possible interactions between genes in the 
network). The long time behavior observed in the network defined on the homogeneous "maximal" graphs 
is featured by the fraction of positive interactions (0 < 77 < 1) allowed between genes. There exists a critical 
value rj c < 1 such that if 77 < r/ c , the oscillations persist in the system, otherwise, when 77 > r/ c , it tends to a 
fixed point (which position in the phase space is determined by the initial conditions and the certain layout 
of switching parameters) . 

In networks defined on the inhomogeneous directed graphs depleted in cycles, no oscillations arise in the 
system even if the negative interactions in between genes present therein in abundance (rj c — 0). For such 
networks, the bidirectional edges (if occur) influence on the dynamics essentially. In particular, if a number 
of edges in the "maximal" graph is bidirectional, oscillations can arise and persist in the system at any low 
rate of negative interactions between genes (rj c — 1). 

Local dynamics observed in the inhomogeneous scalable regulatory networks is less sensitive to the choice 
of initial conditions. The scale free networks demonstrate their high error tolerance. 
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Properties of dynamical networks attract the close attention due to a plenty of their 
practical applications. One class of them is constituted by the growing evolutionary 
networks representing the long time evolution of a genome of living spices or such 
as the World Wide Web where the dynamics and topology evolve synchronously 
according to the external principles of informational safety and economic efficiency 
by the adding of new components at a time. 

Another class is represented by the regulatory networks including a fixed number 
of elements interacting with each other sensitively to the actual position of system 
in its phase space. In the present paper, we study the simplest model of such a 
regulatory network described by a discrete time synchronously updated array of 
coupled structural (topological) and dynamical variables defined at each node of a 
large graph, for both the homogeneous and scalable graphs. 

For small gene expression regulatory networks comprising of just a few elements, 
a direct logical analysis of dynamics is possible, in which their behavior can be 
understood as to be driven by the positive and negative feedback circuits (loops) 
[lj-[2j. However, for large gene expression regulatory networks which can consist of 
thousands of interacting genes, the direct analysis is very difficult because of their 
complexity. Being interested in a qualitative description of dynamical behavior 
exhibited by such large regulatory networks, we turn this problem in a statistical 
way. We suppose that any layout of switching parameters governing the sigmoidal 
type interactions between genes as well as any assignment of gene interaction signs 
(stimulation or inhibition) can be possible with some probability. 

Another important issue discussed in our paper is the influence of global topol- 
ogy onto the local dynamics observed in the large gene expression regulatory net- 
works. Such an effect emerges in many coupled dynamical systems defined on the 
graphs with different topological properties, for instance, in the problem of epi- 
demic spreading [3 J . Being defined on the regular arrays and homogeneous random 
networks, the models of virus spreading predict the existence of a critical spreading 
rate A c > such that the infection spreads and becomes persistent if A > A c and 
dies out fast when A < A c . Recently, it has been shown in [4|-[5j that a variety of 
scale free networks is disposed to the spreading and persistence of infections at 
whatever spreading rate A > the epidemic agents possess that is compatible with 
the data from the experimental epidemiology. 

The dynamical behavior demonstrated by the discrete time coupled map lattice 
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used in the present paper as the models of interacting genes |6J is indeed more 
complicated than that in the probabilistic susceptible - infected - susceptible model 
discussed usually in epidemiology [5]. We show that it is featured by two order 
parameters, that are, the fraction of positive interactions allowed in between genes, 

< r] < 1, and the fraction of bidirectional edges < < 1 presented in the 
"maximal" graph. 

In the large homogeneous regulatory networks, like those defined on the fully 
connected graphs or the regular random graphs, in which all edges are consid- 
ered as bidirectional, the critical fraction of positive interactions in between genes 
at which the oscillations arise and persist in the system is r\ c < 1. In the directed 
inhomogeneous networks, like those defined on the directed scale free graphs, oscil- 
lations die out fast even if the negative interactions between genes present therein 
in abundance (r] c = 0). However, oscillations arise at any low rate of negative in- 
teractions between genes (i] c = 1) provided the "maximal" graph has a number of 
bidirectional edges. Bidirectional edges effectively increase the number of circuits 
presented in the scale free network that is the source of oscillatory behavior. 

The proposed approach could help in understanding the behavior of large gene 
expression regulatory networks for lack of actual empirical data concerning the 
large genomes of living organisms. 

1 Introduction 

An impressive progress in understanding of specific biological regulatory mechanisms which 
play an important role in the way numerous molecular components interact have been made 
recently. This would be a key to control the development and physiology of a whole living 
organism. 

Nevertheless, the behavior of large gene expression regulatory networks is still far from 
being understood mainly because of two reasons: first, up to now, a little is known of the global 
structure of genetic networks. There is still no a common opinion on whether they are organized 
hierarchically (say, as a highly inhomogeneous scalable metabolic network as reported in [HE]) 
or contains a plenty of cross interactions, might be of quite irregular structure, organized in a 
form of connected sets of small sub-networks. Second, the relations between the global structure 
of networks and their local dynamical properties are also still unclear. 

A useful approach to the regulatory networks comprising of just a few elements consists 
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of modelling their interactions by Boolean equations. In this context, feedback circuits (the 
circular sequences of interactions) have been shown to play the key dynamical roles: whereas 
positive circuits are able to generate multistationarity, negative circuits may generate oscillatory 
behavior Genetic networks are represented by the fully connected Boolean networks where 
each element interacts with all elements including itself. A feedback circuit can be formally 
defined as a combination of terms of the Jacobian matrix of the system, with indices forming 
a circular permutation. Flexibility in network design is introduced by the use of Boolean 
parameters, one associated with each interaction of group of interactions affecting a given 
element. Within this formalism, a feedback circuit will generate its typical dynamical behavior 
(either stationary or oscillating) only for appropriate values of some of its logical parameters, 

CHI- 

Most often in biology, the interactions between various molecular components can have a 
definite sign. For any circuit, one can easily check that each element exerts an indirect effect 
on itself which has the same sign for all elements of the circuit, leading to the definition of 
the "circuit sign". In fact, this sign only depends on the parity of the number of negative 
interactions involved in the circuit: if this number is even, then the circuit is positive; if this 
number is odd, then the circuit is negative, 

What makes these concepts important is that specific biological and dynamical properties 
can be associated with each of theses two classes of feedback circuits. The relation between 
the presence of positive feedback loops and the occurrence of multiple states of gene expression 
has been at a focus of investigations for several years (see ^3] and references therein). In 
particular, it has been proven that the presence of positive loop(s) is a necessary condition for 
multistationarity, and the negative circuits (with two or more elements) are required for the 
stable periodicity of behavior, ^U|- Biologically, this means that positive circuits are required 
for different iative decisions and negative circuits are for the homeostasis, [TT]-|13j. 

Nevertheless, for the large regulatory networks comprising of many hundreds or even thou- 
sands of elements, the detailed logical analysis of possible feedback circuits seems to be impos- 
sible now, since the effect that the element can put on itself indirectly, in large regulatory net- 
works, may follow along many different pathways (indeed, if the interaction network is enough 
dense) engaging probably a plenty of distinct feedback structures at once. Furthermore, nu- 
merical observations over the various large discrete time regulatory networks convinced us that 
the " functionable" circuits and the rest "passive" elements are tightly related to each other in 
a way that they play the role of " stabilizers" for the active circuits: a dislocation made to an 
element of the inactive circuit may stampede a change to some "functionable" circuits and even 
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cause they dissolve. 

In the present paper, we focus our attention on the large gene expression regulatory networks 
defined on some "maximal" graph G. To get a qualitative understanding of their dynamical 
behavior, we consider a statistical ensemble of such regulatory networks, in which the thresholds 
assigned to each pairwise gene interaction and its sign are considered as the discrete random 
variables taking their values in accordance to some statistical laws. Let us note that at present 
the values of regulatory parameters driving the behavior of actual gene expression regulatory 
networks are mostly unknown, so that it would be interesting to test the sensitivity of local 
dynamics observed in large regulatory networks to the random change of switching parameters, 
for the different types of such networks. 

Starting from the fixed initial conditions, in large gene expression regulatory networks, we 
shuffled these switching parameters randomly and, otherwise, randomized the initial conditions 
for a fixed layout of thresholds and interaction signs. The long time dynamical behavior ob- 
served in such a statistical ensemble of large regulatory networks depends essentially upon the 
topology of underlying " maximal" graph including all possible pairwise interactions in between 
genes of the given set. Short transient processes arisen in such systems at the onset of simu- 
lations conclude into the statistically stable behavior, that is, either a stationary configuration 
or the multi-periodic oscillations occupying up to a half of system size, in the homogeneous 
regulatory networks having a plenty of negative interactions. 

In contrast to the spreading of chaotic state over the regular and random arrays of piecewise 
linear and logistic discrete time coupled maps studied extensively in the last decade |22.-|24j, 
oscillations arisen in the discrete time regulatory networks do not propagate over the whole 
system and are bounded merely to the oscillating domains. Lack of negative interactions and 
directed cycles in the networks brings it into one of fixed points which position in the phase 
space of system depends upon the certain choice of initial conditions and the layouts of switching 
parameters. The structure of active subgraph of G in the homogeneous regulatory networks 
settled at a fixed point, resembles that one of Erdos-Renyi's random graphs, jEj. 

The plan of paper is following, in Sec. 2, we define the models of large synchronized regulatory 
networks defined on the both homogeneous and inhomogeneous "maximal" graphs where the 
switching parameters are shuffled randomly. In Sec. 3, we present the results of numerical 
simulations on large regulatory networks. In Sec. 4, we introduce the mean field approach 
to the stochastic ensembles of large discrete time regulatory gene expression networks with 
randomly shuffled thresholds of type interactions between genes and their interaction signs. 
Then we conclude in the last section. 
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2 Definition of models 



We define the regulatory gene expression network 1Z on the directed graph G with the set of 
nodes V = {vx, . . . v N } connected by the edges G £ C V x V representing the action of gene 
i onto gene j (the self-action of genes is possible and corresponds to the loops in G). We call 
G as a maximal graph since it contains all possible interactions in between genes of the given 
set. 

The regulatory principle of gene expression networks is that the protein synthesis rate of 
a gene is affected (either stimulated or inhibited) by the proteins synthesized by other genes 
provided their instantaneous concentrations are below (or over) some threshold values. We 
assign the positive sign SV,- = +1 to an interaction if % stimulates the synthesis of protein in 
j, and the negative sign Sij = — 1 otherwise. Indeed, in the real genome, the rate of protein 
synthesis varies from pair to pair of interacting genes, however, for a simplicity, in the present 
paper, we suppose that all interactions between genes are of equal effectiveness, so that all 
edges presented in G have the equal weight. 

The maximal graph G is specified by its adjacency matrix A such that = 1 if G G, 
and A\- = otherwise. Since the effect of interaction between two genes can be negligible at 
time t > (that is, this interaction is switched off at time t), one can define an active subgraph, 
G* C G including all interactions efficient at time t > specified by the instantaneous adjacency 
matrix A*. 

Following [Hj, in the present paper, we consider the synchronized model of gene-gene inter- 
actions, that is, time is discrete and the state of system at time t + 1 is a function of its state at 
time t, in the form of a coupled map lattice. For each gene i, we define two variables: x\ G [0, 1], 
the relative concentration of protein expressing it at time t > 0, and y\ G [0, 1], the exertion of 
gene i, that is, an effective action of other genes onto i at time t which depends upon the their 
relative concentrations x^ at same time t. 

In the homogeneous regulatory networks (like those defined on the complete graphs G(N) 
or on the random regular graphs G(N,r) with the fixed connectivity r > 3, ^3]) the exertion 
y\ can be defined as the fraction of active incoming edges (i.e., the actions of other genes onto 
i) at time t > 0, 



3=1 I 3=1 

The elements of A* are updated synchronously at each time step, in accordance to the current 
values of x* e [0,1]^, 




(1) 




(2) 
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in which the step function 9{z) = 1, for z > 0, and 9(z) = 0, for z < 0, T^- e [0, 1] is the 
threshold value for the action of i onto j, and Sij = ±1 is the sign of interaction. In accordance 
to P|). the interaction « — >■ j is active at time £ if either x\ > Tij for = +1 or < Tjj for 

= — 1. In the case of x\ = Ty, we suppose that the action is active provided that = +1. 

Then, the discrete time synchronous coupling, 



generates the flow 0* in the phase space [0, 1]^ x [0, 1]^ transforming the initial point (x°,y°) 
into (x*, y*) for some £ > 0. The parameter < a < 1 is a positive constant expressing the decay 
rate of protein concentrations, the second term in (J3J) describes the rates of protein synthesis. 
The protein decay rate a determines a time scale in the system, £ — > £(1 — a), and does not 
affect its stable asymptotic behavior. For the first time, such a model has been discussed in jS] 
where the dynamics of several low dimensional examples of genetic expression networks have 
been considered. 

It is important to note that the precise values of switching parameters and are 
at present unknown for the most of pairwise interactions of genes in the genomes of living 
organisms. Being interested in a qualitative behavior exhibiting by the large gene expression 
regulatory networks, we consider these switching parameters as the discrete random variables, 
i.e. we suppose that any layout of them is possible with some nonzero probability, and study 
the statistical behavior of the flow cf) 1 over the ensemble of such networks. Namely, we suppose 
that the threshold assigned to a pairwise interaction of genes can take one of n < N 2 possible 
values, < T\ < . . . < T n < 1, distributed over the unit interval. We also suppose that the 
interaction sign takes the value +1 with some probability < 77 < 1, and = — 1 with the 
probability (1 — 77). In the numerical simulations, we checked out 500 different random strings 
of initial conditions x° for each of 500 random layouts of and S^. 

Dynamical behavior exhibited by the model with the random layouts of switching param- 
eters depends very much on the topology of maximal graph. Many real-world networks and 
gene expression regulatory networks, in particular, can be represented by the inhomogeneous 
graphs of very complex, may be irregular structure having a large number of nodes. They can 
be considered as the quasi-random graphs characterized by some "statistical" properties (for 
instance, by the probability distributions of degrees of their nodes) [T7j . 

It is understood jH] that the formation of feedback circuits driving the dynamical behavior 
in the regulatory networks is closely related to the occurrence of cycles in the underlying "max- 
imal" graphs. In particular, it restricts the application of some random processes (algorithms) 
used for the generation of scale free graphs appropriate for such simulations. For instance, it 




< a < 1, 



(3) 
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is obvious that the popular preferential attachment algorithm proposed in [T^j, in which ver- 
tices are added to the graph one at a time and connected to a fixed number of earlier vertices, 
selected with probabilities proportional to their degrees, is not suitable for the generation of 
such a "maximal" graph, for the gene expression regulatory networks, since it creates a tree 
like structure with no cycles. In the present paper, we use another graph generating algorithm 
(details are given in the Appendix), proposed in and then used in j^j for the study of 
epidemic spreading over the scale free networks. It creates a directed graph with a bounded 
number of cycles for which both the incoming and outgoing degree statistics are scale free. 

In the gene expression regulatory networks defined on the inhomogeneous graphs, the incom- 
ing and outgoing node degrees may differ (like in the directed random scale free graphs), so that 
the definition (JTJ of exertions y\ has to be improved, since, first, some nodes in such graphs may 
have only the outgoing edges (regulating genes), Kj~ = J2f=i ^ij > 0^ while Kf = J2f=i = 0- 
Second, the stability of numerical scheme which ensures that x\ G [0, 1] for any t > requires 
the factor (Kf + Kj~) before the protein synthesis term in 

In the case of directed scale free "maximal" graph, the behavior of the minimally required 
extension of model (JIJ, 

spN At 

y% sr^N 40 , 4O ' ^ ' 

is not very interesting from the dynamical point of view, since all nodes with Y^ =x A\ { < Kf = 
are decoupled from the rest of network, and there are also many nodes (hubs or regulating 
genes) with an excessive number of outgoing edges, K± ^> 1, for which yj <C 1, so that they 
are also effectively decoupled from other genes in the network. 

Given a scale free random graph with the probability degree distributions P(X ± ) ~ (K^)^ 1 
where 7 > 1, then the probability to observe the exertion y l > 1/2 in Q scales with K~ 
like P(y* > 1/2) < (i^-)- 2 ^ 1 YZ=^- l ) m ( m + lY l < h for K ~ > h and decreases rapidly 
with 7. As a result, the most of protein concentrations in the model (jHJ) with the exertion 
defined on the directed scale free graphs are driven by their own decays and got fixed fast close 
to x» = 0. In particulary, one cannot observe oscillations (as a limiting stable behavior) in the 
regulatory gene expression networks (jUEJ) defined on such graphs with 7 > 1, for any layout of 
switching parameters and for any assignment of interaction signs. 

The occurrence of bidirectional edges (when both genes can act onto each other simulta- 
neously) in the highly inhomogeneous scalable graphs can change dramatically the dynamical 
behavior of large regulatory networks defined on them. To quantify the value of bidirectional 
edges in the given "maximal" graph, we introduce the parameter < fi < 1, that is, the fraction 
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of such edges among all edges of G. 

We have performed the numerical simulations for both the homogeneous graphs (the com- 
plete graphs and the "enough dense" regular random graphs comprising of 10 3 nodes) assuming 
all edges in them to be bidirectional, with no loops, and the highly inhomogeneous scalable 
graphs, the directed scale free graphs G(10 3 , 2.2) such that both probabilities that a node has 
precisely K + incoming edges and K_ outgoing edges follow the power law, P(K±) ~ K± 2 ' 2 . 
This graph has been reported in [7] as being typical for the metabolic networks of many living 
organisms. In the latter case, we have varied the fraction of bidirectional edges in the whole 
interval < // < 1. 

3 Numerical results 

We monitor the system approaching a statistically stable regime by tracking out its "veloci- 
ties" in the phase space, v* = x m — x* (the rate of protein synthesis) and u l = l^-ij 1 ~ 
/N 2 (the rate of gene exertions) counting the number of edges triggered between G* and G\ G* 
at time t. While studying the model of large regulatory networks (j2HHJ) with the random layouts 
of switching parameters defined on both homogeneous and inhomogeneous graphs consisting of 
10 3 nodes, we have varied the number n of distinct threshold values, < Ti < . . . < T n < 1, 
from several tens to several hundreds changing by this way the coarse graining of phase space. 
We choose the threshold values {Ti, . . . , T n } uniformly distributed (u.d.) over the interval [0, 1]. 
In fact, in the actual simulations, we have assumed the "micro canonical" distribution of thresh- 
olds characterized by the probability density function g (x) = n' 1 X}fc=i S(x — k/n), where 5(x) 
is the Dirac function, so that each "event" = k/n had the same statistical weight, n~ l . 

To display the dependence of dynamical behavior on the random initial conditions or the 
random layouts of switching parameters or both, we have presented the collected statistical data 
either in the form of density plots or the box-plots where a central line of each box shows the 
median (i.e., the middle value collected over 500 random initial conditions or random layouts), 
a lower line of a box shows the first quartile, and an upper line of a box shows the third quartile. 
Half of the difference between the third quartile and the first quartile (the semi-interquartile 
range or the quartile deviation) is a measure of the dispersion of the data. Two lines extending 
from the central box of maximal length 3/2 the interquartile range. 

We commence the analysis with the networks defined on the undirected homogeneous graphs: 
the complete graph and the regular random graph in which the connectivity of each node is 
fixed at some K > 3 to ensure the existence of many Hamilton cycles traversing all nodes in the 
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graph. All nodes of such graphs have the equal connectivities K and each node is connected to 
any other with some probability p > (for the fully connected graph, K = N — 1 and p — 1). 
The behavior of the model (001) defined on them looks essentially similar: after a considerably 
short transient process, the system settles into a statistically stable dynamical regime which 
depends very much upon the fraction of negative interactions allowed in between genes and 
quantified by 1 — rj. 

Provided the system has just a few negative interactions in the network (rj — > 1), it ap- 
proaches one of the fixed points exponentially fast, for any random initial condition x° (see 
Fig. [TJa). Herewith, the rate of synthesis of proteins (|v| 4 ) shows a crossover between two 
consequent phases of transient regime (see Fig. [I] a): before and after the moment when the 
structure of active subgraph G* gets fixed. The exponential decay rate of transients is indepen- 
dent on neither the choice of initial conditions nor the certain layout of switching parameters, 
but the stationary asymptotic configuration x* depends upon both. 

For a given random layout of switching parameters, many different stable fixed points can 
be achieved by the system starting from the different initial conditions x° that is an evidence of 
multistationarity in the model. To get an image of variety of the possible asymptotic stationary 
states, we have counted the occupancy numbers pk = {§i : T k -\ < x^ < T k } of the consequent 
intervals of phase space, A k = [T fc _x,T fc ], formed by 100 threshold values T k u.d. over the 
unit interval (see Fig. QJb). The boxes shown in Fig. ^b present the variations of occupancy 
numbers p k for 500 different random initial conditions x° with 5% of negative interactions 
allowed between genes (rj = 0.95), for some fixed layout of switching parameters. 

The density of possible stationary asymptotic configuration x* depends upon the certain 
layouts of switching parameters for any certain initial string x° that is presented on Fig. |21 A 
patchy structure of graph in Fig.|21a. manifests the multistationarity in the system, meanwhile 
the "clusters" formed by the merged patches show that the limiting stationary configurations x* 
are sensitive to the layout of switching parameters and can move gradually as these parameters 
shuffle. 

Shuffling of switching parameters mixes up the orbits of deterministic system intensively, 
as a result each x\ spreads out fast with time over the whole unit interval. In Fig. |21b, we 
have sketched the density plot of possible values of x\ (for % = 77 in G(10 3 )) vs. time in 30 
consequent time steps (long enough to achieve a fixed point) starting from the initial value 
x 7 = 0.175. 

It is worth to mention that at any fixed point the active subgraph G* C G constitutes a 
random graph ("half-dense" in comparison with G). In Fig. |31 we have shown the probability 
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degree distributions (the circles are for the incoming degrees K + , and the diamonds are for the 
outgoing degrees K_) for the nodes of active subgraph formed at a fixed point of the model 
defined on the fully connected graph G(10 3 ) with 1% of negative interactions allowed 
between genes. The solid line on Fig. |3] displays the Gaussian probability degree distribution 
which is typical for the Erdos and Renyi random graphs, |14j . 

One can say that in the homogeneous regulatory networks when the positive interaction 
between genes prevail in the system, its dynamical behavior is dominated by the positive feedback 
circuits responsible for a number of asymptotically stable states (fixed points). Thereat, the 
strain of negative feedback circuits related to just a few negative interactions is statistically 
negligible. 

With increasing percentage of negative interactions allowed in the system up to approxi- 
mately 10%, it exhibits a complicated spatiotemporal behavior where the domains of genes 
with the stationary concentrations of proteins coexist and interleave with those of periodically 
oscillating concentrations. In contrast to the spatiotemporal intermittency observed in the syn- 
chronously updated discrete time extended dynamical systems defined on the various regular 
arrays [221 E3] and on the regular random graphs [21] , the dynamical state (oscillations) does 
not propagate with time throughout the regulatory network. 

The oscillating domains arisen in the homogeneous regulatory networks are bounded by 
the genes whose oscillation amplitudes are insufficient for their protein concentrations to cross 
the next thresholds. In the large enough homogeneous regulatory networks, a turnover of 
nodes engaged into the oscillating domains happens occasionally. The averaged number of 
nodes joined the oscillating domains at a time increases as the fraction of negative interactions 
allowed in the model grows up. In Fig. HJa, we have displayed the decreasing and vanishing of 
fractions N osc /N of oscillating nodes with rj, in the model defined on the fully connected graph 
(with bidirectional edges). Boxes show the fluctuations of these fractions in the ensemble of 500 
different random layouts of switching parameters and random initial conditions. Bold points 
stand for the means of collected data. 

The solid line in Fig. 0]a is the Gaussian curve, N osc / N ~ exp(— rf /2a)/ V2 ir a 2 , with 
a ~ 0.555 fitting the data well. The direct logical analysis of low dimensional regulatory 
networks (see, for example ^H]) relates oscillations to the dynamical patterns generated by 
the "functional" negative feedback circuits. A feedback circuit exhibits its typical dynamical 
behavior only for appropriate values of some of the logical parameters. The graph sketched on 
Fig. EJa demonstrates that the probability to observe a "functional" negative feedback circuit, 
in large homogeneous regulatory networks with the randomly shuffled switching parameters, is 
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close to the normal statistical law with regard to rj. 

Fig. EJb displays the distribution of nodes changing their protein concentrations periodically 
vs. the periods of such changes observed in the large homogeneous regulatory networks. The 
distribution in Fig. |3]b counts all such nodes disregarding for the amplitudes of changes. Each 
node has been counted in the distribution just once under the minimal period of oscillations 
it exhibits. As usual, boxes represent the fluctuations of numbers of oscillating genes N(t) 
over the ensemble of 500 different random layouts of switching parameters and random initial 
conditions. The distribution has a maximum at r = 3 independently upon the initial conditions 
and the layouts of switching parameters. 

Formation of oscillating domains for 77 1 (when the negative interactions present in abun- 
dance) comes along with the synchronization of the rest of system at x = 1/2 (see the profile 
for the occupancy number pk in Fig. EJa). This synchronization looks essentially insensitive to 
the initial conditions (the boxes on the graph are almost imperceptible) . It gives an impression 
that when the dynamical behavior is obviously driven by the negative feedback circuits, and 
oscillations of protein concentrations can occupy up to a half of nodes in the network, just a 
few of them actually change the instantaneous structure of active subgraph G*. 

On Fig. Ob, we have shown the behavior of the phase space velocities, log |v*| and log \u l \, 
vs. logt characterizing the transient processes in the system for the model defined on the 
fully connected graph G(10 3 ), with 100 distinct thresholds u.d. over the unit interval, with 
rj = 0.5 (1/2 of interactions are negative). When the negative interactions prevail in the system, 
these velocities decay much slower than the exponentially fast transients shown in Fig. ^a and 
asymptotically tend to a power law as rj — ■> 0. Moreover, they do not extinguish eventually and 
bring in oscillations in short time. 

Stable patterns of statistical behavior are insensitive to the random layouts of thresholds 
and assignments of interaction signs. Starting from some fixed initial string x°, we have shuffled 
the switching parameters in the model defined on G(10 3 ) (with bidirectional edges) for rj = 0.5 
and displayed the data for x\ (after the stable behavior had been achieved) for the first 150 
nodes (see Fig. El a). The time evolution of these density distributions indicating oscillations 
is illustrated by the beatings shown in Fig. Elb in 30 consequent time steps taken over the 
ensemble of 500 different random layouts of switching parameters. 

It is important to note that for any value of r\ it is always less than 1/2 of nodes (achieved 
only if all interactions between genes are negative) that are engaged into oscillations in the 
statistically stable regime. Herewith, the oscillating protein concentrations x\ for the most of 
nodes are still bounded within their intervals = [T^-i, Tk] and do not change the occupancy 
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numbers pk- Direct investigation of active subgraphs in the stable oscillating regimes had shown 
that for any random initial condition the notable oscillations of protein concentrations resulting 
in changes to active subgraph are confined to less than 3% of nodes. 

To illustrate the multi-periodicity and regularity of oscillations arisen in the model (J2H3J) 
defined on G(10 3 ), we have presented the return maps expressing log|v t+1 | vs. log | v* | and 
log |w t+1 | vs. log for 77 = 0.5 in Fig. [7|a and for 77 = 0.1 in Fig. [7|b. Detailed observations 
of the fine structure of return maps in Fig. 0a show that when the system has a parity be- 
tween negative and positive interactions, it exhibits oscillations of shortest periods, r = 2,3,4. 
Oscillations of many different periods arise if the negative interactions prevail in the system. 

In a conclusion, when the fraction of positive interactions between genes in the regulatory 
network defined on the homogeneous graphs decreases down to the critical value r] c ~ 0.92, 
the negative feedback circuits influence the dynamical behavior substantially that is revealed 
in persistent oscillations arisen in the system at the end of transient processes. 

We have also studied the behavior of model (j2HHJ) defined on the inhomogeneous scale free 
graph G(10 3 , 2.2) with the exertion (@J), in which both statistics of incoming and outgoing edges 
follow the power law P(i ; r ± ) ~ (_fr ± ) — 2-2 shown on Fig. [HI It has been generated in accordance 
to the algorithm given in Appendix A. 

The behavior demonstrated by the dynamical systems defined on scale free graphs depends 
essentially on their detailed topology and can be dramatically different for the graphs gener- 
ated in accordance to distinct algorithms even though they enjoy the same probability degree 
statistics Individual structural properties of scale free graphs is characterized by the certain 
pair-formation process, in which each vertex v of degree K chooses a set of partners according 
to a specified K— dependent rule describing its "preferential choice". Directed scale free graphs 
have usually a few cycles (since cycles imply a statistical parity between numbers of incoming 
and outgoing edges and therefore cannot play a key role in the structures of quasi-random 
scalable graphs). 

The feedback circuits formation, in scale free graphs, relays primarily upon the bidirectional 
edges connecting pairs of mutually interacting genes. In Fig. we have presented a phase 
diagram displaying the appearance of persistent oscillations in the model (j2HHJ) defined on the 
inhomogeneous scale free graph G(10 3 , 2.2) for different fractions of negative interactions rj 
and bidirectional edges /i. It shows that for the enough fractions of bidirectional edges, in 
the scalable regulatory network, oscillations can persist even if there is just a few negative 
interactions allowed in between genes of the network (j] c = 1). 

We have studied the statistic of oscillations, in the limiting case of /1 = 1, when the scale free 
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graph can be considered as undirected. Fig. [TUla displays the decrease of oscillating domains 
(with a rate close to linear) in the scalable regulatory network. Fluctuations of data collected 
over the ensemble of 500 different random layouts of switching parameters and initial conditions 
at a = 0.74 are shown by boxes. Let us note that Fig. EBa reveals the similar effect of 
global scalable topology of network onto the local dynamics of nodes that has been reported 
recently in on the problem of virus spreading in the undirected scalable networks. In 

the homogeneous networks, the spreading of a dynamical state is a threshold phenomenon: 
it occupies a valuable fraction of network nodes as the control parameter determining the 
formation rate of the state exceeds some critical value. Alternatively, in the scalable networks, 
the dynamical state spreads and persists at whatever value of control parameter. 

Fig. EBb presents the distributions of nodes with periodically changing protein concentra- 
tions vs. the periods of changes in the undirected scale free graph G(10 3 ,2.2) for rj = 0.1 (the 
upper profile) and for 77 = 0.5 (the lower profile) taken at a = 0.74, // = 1. Distributions 
count all nodes of network which change their protein concentrations periodically disregarding 
for the amplitudes of these oscillations. Each node has been counted in the distribution just 
once under the minimal period of oscillations it exhibits. It is seen that the maximal number 
of nodes demonstrates the minimal period of changes r = 2 that is not a surprise since the 
most of negative feedback loops of minimal periods in the scalable networks relays upon the 
bidirectional edges. 

The statistics of stable long time behavior observed in the scalable regulatory networks 
can be characterized by the occupancy numbers pt of the consequent intervals of phase space, 
Afc = [Tfc_i,Tfc]. In Fig. El we have shown the occupancy numbers of the model defined on 
the undirected scale free graph G(10 3 , 2.2) with a given configuration of 50 thresholds u.d. 
over the unit interval for two opposite values of rj. Fluctuations shown by the boxes reveal the 
dependence of the occupancy number upon the certain choice of initial conditions and layouts 
of switching parameters. It is important to note the high error tolerance of scalable regulatory 
network: in its phase space, the intervals exist which stay void (p& = 0, even though some 
nodes had initially been settled into these intervals) and for which pk is got fixed independently 
upon the layout of switching parameters and initial conditions. 

When the negative interactions present in abundance (77 <C 1), the valuable fraction of nodes, 
in the scalable regulatory network, synchronizes either in the first (next to x — 0) or in the 
last (next to x — 1) intervals of phase space. The nodes demonstrating oscillations of protein 
concentrations (about N/2 total as 77 — ► 0) are scattered in the middle range of phase space. 



14 



4 Mean field approach to the large regulatory networks 



The aim of present section is to introduce a "mean field" approach to the large regulatory 
networks with randomly shuffled switching parameters. 

In accordance to (J2H21), the j-th protein synthesis rate depends upon the exertions t/i of all 
genes acting on it, that is, the fraction of active arcs incident to Vj at time t. While shuffling 
randomly the switching parameters in the large regulatory network, we suppose that the positive 
sign S = +1 is selected for the action % — ► j with the probability < r\ < 1, and the interaction 
threshold value equals to some T k (k = 1, . . . , n) chosen with some probability < Pq. < 1, 
such that Y^k=\ Pik — 1 f° r eacri gene i acting on j. 

The distribution of values T k over the unit interval can be defined by the set of integrable 
functions g k ■ [0, 1] — > 1R + satisfying the normalization condition g k (z) dz = 1 and such that 
the probability that the threshold value T k chosen for the interaction does not exceed x reads 
as 

P(T fc < x) = I g k (z) dz. 



Then the probability that an outgoing arc is active at time t equals to 



A v( x l) = J2 Pik y J 1 9k(z) dz + (l- V ) J 9k(z) dz) , (5) 



k=l 

so that the expected exertion of the gene j at time t is 

N IN 

^) = E4 A ^)/E4- ( 6 ) 

i=l / i=l 

The main idea of mean field approach to the large regulatory networks with shuffled switching 
parameters is to replace the exact exertion values y\ in (J3J) with their expectations (JUJ), 

(x* +1 ) = a(x t ) + (l-a)(y*), (7) 

in which the expected exertion (y*) is calculated in accordance to (JHJ). 

Generally speaking, the problem (J7|) is as complicated as the original one. However, in some 
cases it can be more "illustrative" revealing the entire dynamical mechanisms of the system. 
In the simplest case of uniform probabilities = 1/n, the system of N coupled equations ((7j) 
is linearized for the uniform distributions, — 1. After the exponentially fast extinguishing 
of transient processes, the system (J7|) with constant P^ and g k is settled into the uniformly 
stationary configuration, 

rj)n 



X*i 



(n + l-2rj)(N- I)' 
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independently upon the initial conditions. 

For the "microcanonical" distributions of thresholds, = 5(x — T^), with the randomly 

chosen threshold values Tk u.d. over the unit interval and Pik = 1/n, the dynamics of expecta- 
tions over the ensemble of large networks (N ^> 1) with randomly shuffled switching parameters 
is given by 

= a(xf) + (1 - a) [(1 - V )f(T m ) + V (l- f(T m ))] (8) 

where the "mean field" /*(T m ) = N(x t > T m )/ N is the instantaneous fraction of genes with 
concentrations x l > T m . The "mean" threshold value T m is determined upon the certain set 
of randomly chosen thresholds {T^} by the equation (|SJ) in the limit N ^> 1. In this case, 
the equation (jHJ) represents a system with a feedback loop. Herewith, depending upon the 
value of rj and T m it can be either positive or negative. The system governed by (JSJ) is getting 
synchronized very fast, in particular, it exhibits the synchronous oscillations of different periods 
if rj < T m < 1 — T), for r] < 1/2, and T m < 1 — r], T m > rj, for rj > 1/2. 

In the case of inhomogeneous sets Pa-, for N ^> 1, the equation (jHJ) can be generalized to 

(x? 1 ) = a(xf) + (1 - a) [(1 - V )f(T mi ) + V (l- ?{T mi ))] (9) 

where T mi is the local "mean" threshold related to gene i. The next step can be done if one 
supposes some symmetry properties for the sets of Pik that dramatically reduces the number of 
independent local "mean" thresholds T mi to just a few. In such a case, the entire system can 
be considered as a set of several coupled positive and negative feedback loops whose dynamical 
behavior can be a subject of detailed analysis. In particular, in the system with several cou- 
pled negative feedback circuits characterized by the distinct thresholds T mi , the synchronous 
oscillations may persist for rj > 0. 

5 Conclusion 

In the present paper, we have studied the behavior of large discrete time regulatory networks 
defined on the homogeneous and scalable maximal graphs. The "traditional" approach to 
the complex regulatory networks comprising of a few elements reduces their modelling to the 
Boolean equations. In this context, the concept of "feedback circuits" has been discussed 
extensively. However, for the large regulatory networks, a direct logical analysis seems now 
impossible because of their complexity. 

To get a qualitative insight into the behavior of regulatory networks of thousands elements, 
we study how sensitive it is upon the shuffling of switching parameters, that are, the interaction 
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signs Sij = ±1 and the interaction threshold values Ty. Starting from a fixed initial condi- 
tions, we have shuffled these switching parameters randomly and then randomized the initial 
conditions for a given layout of thresholds and interaction signs in the network. 

The dynamical behavior observed in such a system crucially depends upon the topology of 
maximal graph including all possible interactions of genes in the given network. 

We have found that, in the homogeneous regulatory networks, the parameter < 77 < 1 bal- 
ancing the positive and negative interactions allowed between genes is a convenient parameter 
featuring their dynamical behavior. 

When the positive interaction between genes prevail in the system, its dynamical behavior is 
dominated by the positive feedback circuits responsible for a number of asymptotically stable 
states, at the same time the influence of negative feedback circuits onto the dynamical behavior 
of system is negligible. The structure of active subgraph in the regulatory network settled at 
a fixed point is close to the random graph of Erdos and Renyi, ^3]. With the increase of 
percentage of negative interactions up to approximately 10%, the system exhibits a complicated 
spatiotemporal behavior where the domains of genes with the stationary concentrations of 
proteins coexist and interleave with those of periodically oscillating concentrations. 

The feedback circuits formation, in the regulatory networks defined on the scale free graphs, 
relays primarily upon the bidirectional edges connecting pairs of mutually interacting genes. 
For the enough fractions of bidirectional edges, in the scalable regulatory networks, oscillations 
can persist even if there is just a few negative interactions allowed between genes of the network. 

The model of gene expression regulatory networks which we have considered reveals the 
similar effect of global scale free topology on the local dynamics of nodes that has been reported 
recently in on the fractions of infected agents in the undirected scalable networks vs. the 

effective spreading rates A. In the homogeneous networks, a dynamical state occupies a valuable 
fraction of network nodes as a control parameter determining the formation rate of the state 
exceeds some critical value. In contrast, in the scalable networks, the dynamical state spreads 
and persists at whatever value of control parameter. 

In the phase space of scalable regulatory networks, the intervals exist which stay void (the 
occupancy numbers pk = 0) and for which p^ is got fixed independently upon the layout of 
switching parameters and initial conditions. As usual, the scale free networks demonstrate their 
high error tolerance. 

The important unsolved problem deserving a thorough investigation is that of dynamics of 
oscillating domains in the graph. Due to the presence of translational invariance (in average) 
in the homogeneous large regulatory gene expression networks, the oscillating domains can 
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probably move in the graph by the interchange of nodes with the stationary part of the network. 
We can also suggest that perhaps, in the inhomogeneous large regulatory gene expression 
networks defined on the scale free graphs, the motion of oscillating domains proceeds not 
in the "real space" of the graph, but rather in between different classes of nodes having the 
same connectivities. For instance, we have beheld that while propagating throughout a scale 
free graph, the dynamical state occupies at first the nodes with minimal connectivities, then 
those of maximal connectivities (hubs) and eventually the nodes of intermediate connectivities. 
Similar dynamical properties of the scale free networks has been reported recently in |5J. We 
shall discuss this issue in the forthcoming publications. 
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7 Appendix A: 

The large scale metabolic networks have been studied extensively in It has been shown 
that for many different living organisms the probability degree distributions for both incoming 
and outgoing edges of a node in the interaction graphs follow the power law, P{K) oc K~ 2 ' 2 , 
where K is the connectivity of a node. A flexible algorithm generating scale free graphs based on 
the principle of evolutionary selection of a common large-scale structure of biological networks 
incompatible with the preferential attachment approach ( see jS]), has been discussed in 
PU] recently. A random procedure which we use to generate the scale free graph draws back to 
the "toy" model for a system being at a threshold of stability reported in [21]. Here we briefly 
explain this algorithm for convenience of the readers. 

One consider three random variables x, y, and z that are the real numbers distributed in 
accordance to the distributions /, g, and v within the unit interval [0,1]. We assume that x 
represents the current performance of a biological network (say, the protein-protein interaction 
map), while y and z are the thresholds for outgoing and incoming edges respectively. The net- 
work is supposed to be stable until x < y and x < z, and is condemned otherwise. Fluctuations 
of thresholds reflect the changes of an environment. 
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The random process begins on the set of N vertices with no edges at time 0, at a chosen 
vertex i. Given two fixed numbers fi G [0, 1] and v e [0, 1], the variable x is chosen with respect 
to pdf /, y is chosen with pdf g, and z is chosen with pdf v , we draw e^- edge outgoing from 
% vertex and entering j vertex if x < y and x < z and continue the process to time t = 1. 
Otherwise, if x > y (x > z), the process moves to other vertices having no outgoing (incoming) 
links yet. 

At time t > 1, one of the three events happens: 

i) with probability fi, the random variable x is chosen with pdf / but the thresholds y and 
z keep their values they had at time t — 1. 

ii) with probability 1 — /i, the random variable x is chosen with pdf /, and the thresholds y 
and z are chosen with pdf g and v respectively. 

iii) with probability u, the random variable x is chosen with pdf /, and the threshold z is 
chosen with pdf v but the threshold y keeps the value it had at time t — 1. 

If x > y, the process stops at i vertex and then starts at some other vertex having no 
outgoing edges yet. If a; > z, the accepting vertex j is blocked and does not admit any more 
incoming link (provided it has any). If x < y and x < z, the process continues at the same 
vertex i and goes to time t + 1. 

It has been shown in [20] that the above model exhibits a multi-variant behavior depending 
on the probability distribution functions /, g, and v chosen and values of relative frequencies 
\i and v. In particular, if v — 0, both thresholds y and z have synchronized dynamics, and 
sliding the value of fi form to 1, one can tune the statistics of out-degrees and in-degrees 
simultaneously out from the pure exponential decay (for fi = 0) to the power laws (at /j, — 1 ) 
provided /, g, and v belong to the class of power law functions. For instance, by choosing the 
probability distribution functions in the following forms 



f(u) = (1 + a)u a , a>-l, 
V (u) = g(u) = (l + f3)(l-uy, (3>-l, 



(10) 



one obtains that 



— n±ms^±^ (i + o (i)) . pi, 

For different values of (3, the exponent of the threshold distribution, one gets all possible 
power law decays of p^ = i(K). Notice that the exponent 7 = 2 + f3 characterizing the decay 
of Pyt=\(K) is independent of the distribution f(u) of the state variable x. In the uncorrelated 
case, /j, = 0, the degree distribution functions decays exponentially (for instance, = 2~ K 
for / = g = v = 1) [20] • For the intermediate values of /1, the decay rate is mixed. 
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The preference attachment matrix Hsk which elements are the probabilities that a vertex 
with degree S is connected to a vertex having degree K, for a scale free graph, generated in 
accordance with the above algorithm depend only of one variable K, jj], 



Expanding the binomial in the above equation, one gets the leading term oc (K/N — I) 13 . 
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8 Figures 
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Figure 1 : The behavior of model (121,'i^ for a given random layout of switching parameters defined on the fully 
connected graph G(10 3 ) in the vicinities of its fixed points, a). The exponentially fast approach of the rate of 
synthesis |v*| = |x* +1 — x*| and u* = — \ to zero for the model with a fixed configuration of 

50 thresholds uniformly distributed (u.d.) over the unit interval and the fixed assignment of S = —1 to 1% of 
interactions, 77 = 0.99, the protein decay rate a = 0.7. b). The occupancy number pk of the consequent intervals 
Afc = [Tfc_i,Tfe] at the different fixed points of model with a given configuration of 100 thresholds u.d. over 
the unit interval, r\ = 0.95 (5% of negative interactions), a = 0.7. The strong fluctuations shown by the boxes 
reveal the dependence of the occupancy numbers upon the certain choice of initial conditions; data has been 
collected over 500 random initial strings x°, for a fixed layout of switching parameters. 
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Figure 2: The density of stationary configurations x* at the fixed points depends strongly upon the certain 
layouts of switching parameters, a). The density plot of empirical distributions of stationary values Xi* for 100 
nodes chosen at random among total 10 3 . Stationary configurations have been achieved by the system starting 
from the given string of initial conditions. 50 distinct threshold values u.d. over the unit interval have been 
shuffled (500 different layouts) randomly together with the random assignments of 1% of negative interactions 
(rj = 0.99), a = 0.6. The patchy structure of graph reveals the multistationarity in the system. Some patches 
merge manifesting the sensitivity of stationary configurations to the certain layout of switching parameters, b). 
The density plot of empirical distributions of values x\ 7 vs. time in 30 consequent time steps (long enough to 
achieve a fixed point) starting from x° 7 = 0.175 on the fully connected graph of 10 3 nodes. The random shuffle 
of switching parameters mixes up the orbits of deterministic system \'1W,'<\ . 
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Figure 3: The probability degree distributions of nodes in the active subgraph formed at a fixed point of the 
model (121: 11) defined on the fully connected graph G(10 3 ) with 1% of negative interactions allowed between genes 
(the circles are for the incoming degrees and the diamonds are for the outgoing degrees of nodes). The solid line 
is for the Gaussian probability degree distribution which is typical for the random graphs of Erdos and Renyi, 

m 
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Figure 4: a). Decreasing and vanishing of oscillating domains as r\ approaches rj c ~ 0.92 in the regulatory 
network defined on the complete graph G(10 3 ). Boxes present the fluctuation of data collected over the ensemble 
of 500 different random layouts of switching parameters and initial conditions, a — 0.74. The bold points stand 
for the means. The data are fitted well with the Gaussian curve (the solid line), N osc / N ~ exp(— rj 2 
, in which a ~ 0.555. b). The distribution of nodes with the oscillating protein concentrations over the periods 
of oscillations in the complete graph <G(10 3 ) at rj = 0.1, a = 0.74. Boxes present the fluctuation of data collected 
over the ensemble of 500 different random layouts of switching parameters and initial conditions. 
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Figure 5: a). The occupancy number pk of the consequent intervals = [Tk-i,Tk] is insensitive to the 
initial conditions for the random layout of switching parameters as J) C 1 (computed for the model defined 
on the complete graph G(10 3 ) with 100 distinct thresholds u.d. over the unit interval, 77 = 0.5, a — 0.6). 
When oscillations arise, the values x\ for the most of nodes (which do not engaged into these oscillations) are 
synchronized at 1/2. b.) When the fraction of negative interactions increases, the transient processes (which 
decay following a power law as 77 — ► 0) in the model (ll'loli defined on the homogeneous graphs conclude in 
oscillations. 
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Figure 6: Oscillations arise in the system together with synchronization of the rest of network at 1/2 is 
insensitive to the shuffle of thresholds and interaction signs, a.) The density plot of empirical distributions of 
values x\ for the first 150 nodes of the complete graph G(10 3 ) in the stable oscillating regime, 77 = 0.5, a = 0.6, 
taken over 500 distinct shuffles of switching parameters, b.) The density plot of empirical distributions of 
values £77 vs. time in 30 consequent time steps (t = 100 . . . 130) taken over the ensemble of 500 different 
random layouts of 50 distinct thresholds. 
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Figure 7: The return maps for the model Q2I3JI defined on <G(10 3 ) for the fixed layout of 100 distinct thresholds 
u.d. over the unit interval and fixed random initial conditions x , a — 0.6, reveal the multi-periodicity of 
oscillations risen in the system, a.) for 77 = 0.5, b.) for r\ = 0.1. 
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Figure 8: The scale free graph G(— 2.2, 10 3 ) generated in accordance to the algorithm explained in the Appendix 
A and used for the numerical simulations on the model 1)213(1 . It is characterized by the binomial preference 
attachment Tl{K) ~ (1 — K / N — l) ' 2 . Both probability degree distributions for in-degrees and out-degrees 
follow the power law P(K) oc K~ 2/2 that is reported to be typical for the large scale metabolic networks of 43 
different organisms studied in [Jj. On the figure, we have drawn the vertices with higher connectivities closer to 
the center of the graph and those of lower connectivities are located on the periphery, so that the nodes equally 
distanced from the center have equal connectivities. 
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Figure 9: The phase diagram (77, n) displays the appearance of persistent oscillations in the model ()2l.'{t defined 
on the inhomogeneous scale free graph G(10 3 , 2.2) taken at a = 0.3. Contours bound the regions with the equal 
numbers of oscillating nodes. For small fractions of bidirectional lines fi <C 1, oscillations fade out, while they 
arise when [i > 0.2 persisting at any rate r\ > 0. The details of the phase diagram strongly depends upon the 
precise structure of random scale free graph. 
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Figure 10: a). Decreasing of oscillating domains in the regulatory network defined on the undirected scale free 
graph G(10 3 ,2.2) at /x = 1. Boxes present the fluctuations of data collected over the ensemble of 500 different 
random layouts of switching parameters and initial conditions, a — 0.74. The bold points stand for the means, 
b). The distributions of nodes with oscillating protein concentrations vs. the periods of oscillations in the 
undirected scale free graph G(10 3 , 2.2) at ij = 0.1 (the upper profile) and at r\ — 0.5 (the lower profile) taken at 
a = 0.74 and \i = 1. 



Figure 11: The occupancy numbers pk of the consequent intervals of phase space, Afe = [7fc_i, jy , in the model 
defined on the undirected scale free graph G(10 3 , 2.2) with a given configuration of 50 thresholds u.d. over the 
unit interval at a = 0.7. Fluctuations shown by the boxes reveal the dependence of occupancy numbers upon 
the certain choice of initial conditions and layouts of switching parameters (500 different configurations had 
been used to collect data), a), for i] = 0.99 (1% of negative interactions), b). for rj = 0.01 (99% of negative 
interactions). 
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